home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / raytrace / radiance / simplerd.lha / simplerad / FinalFTP / Light / raypoly.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-05-21  |  10.7 KB  |  347 lines
  1. /**********************************************************************/
  2. /* ray_polygon.c                                                      */
  3. /*                                                                    */
  4. /* Ray / polygon intersection routines                                */
  5. /*                                                                    */
  6. /* Copyright (C) 1992, Bernard Kwok                                   */
  7. /* All rights reserved.                                               */
  8. /* Revision 1.0                                                       */
  9. /* May, 1992                                                          */
  10. /**********************************************************************/
  11. #include <stdio.h>
  12. #include <math.h>
  13. #include "geo.h"
  14. #include "misc.h"
  15. #include "struct.h"
  16. #include "ray.h"
  17. #include "io.h"
  18.  
  19. extern OptionType Option;
  20. extern RayStats_Type RayStats;
  21. extern Polygon *shootPatch;
  22. extern Polygon *recPatch;
  23.  
  24. /**********************************************************************/
  25. #define Positive_side 1
  26. #define Negative_side 0
  27. #define get_side(y) (y>=0.0 ? 1 : 0)         /* Find out side of poly */
  28.  
  29. /**********************************************************************/
  30. /* Find maximal component of a vector */
  31. /**********************************************************************/
  32. #define get_max_component(v) \
  33.   (((fabs((v).x)>fabs((v).y))&&(fabs((v).x)>fabs((v).z)))?0 \
  34.    :(((fabs((v).y)>fabs((v).x))&&(fabs((v).y)>fabs((v).z)))?1:2))
  35.  
  36. /**********************************************************************/
  37. /* Get u and v directions based on 2 points, and normal direction     */
  38. /**********************************************************************/
  39. void get_uv(w,p1,p2,u,v) 
  40.      int w;
  41.      Vector p1, p2;
  42.      double *u, *v;
  43. {
  44.   switch(w){
  45.   case 0 : 
  46.     *u = p1.y-p2.y; *v = p1.z-p2.z; 
  47.     break;
  48.   case 1 : 
  49.     *u = p1.z-p2.z; *v = p1.x-p2.x; 
  50.     break;
  51.   case 2 : 
  52.     *u = p1.x-p2.x; *v = p1.y-p2.y; 
  53.     break;
  54.   default:
  55.     break;
  56.   }
  57. }
  58.  
  59. /**********************************************************************/
  60. /* Test if ray-plane intersection is inside a polygon                 */
  61. /* Project the polygon points on a co-ordinate plane. The projection  */
  62. /* does not preserve area but topology is preserved. Use projection to*/
  63. /* count polygon edge crossings. If the number is odd then the point  */
  64. /* is inside else it is outside.                                      */
  65. /**********************************************************************/
  66. int Inside_Polygon(p,pptr)
  67.      Vector *p;                    /* Intersection point */
  68.      Polygon *pptr;                /* Polygon to test */
  69. {
  70.   int i,ii,max_comp;
  71.   int first_x_status, first_y_status, second_x_status, second_y_status;
  72.   int num_of_crossings;
  73.   double u1, u2, v1, v2;
  74.   int numVert;
  75.  
  76.   /* Find the the component with largest magnitude 
  77.      in the normal to the plane of the polygon. 
  78.      Suppose x component has largest magnitude, then plane is yz
  79.      Compute u, v directions on this plane.
  80.      */
  81.   num_of_crossings = 0;
  82.   numVert = pptr->numVert;  
  83.   max_comp = get_max_component(pptr->normal[0]); 
  84.   get_uv(max_comp,pptr->vert[0]->pos,*p,&u1, &v1);
  85.  
  86.   /* Traverse poly and count crossings */
  87.   first_x_status = (u1 > 0.0);
  88.   first_y_status = (v1 >= 0.0);
  89.  
  90.   for(i=1;i<=numVert;i++) {    
  91.     ii = ((i==numVert) ? 0 : i);
  92.     get_uv(max_comp, pptr->vert[ii]->pos, 
  93.        (*p), &u2, &v2);
  94.     second_x_status = (u2 > 0.0);
  95.     second_y_status = (v2 >= 0.0);
  96.  
  97.     /* Count crossings */
  98.     if (first_y_status != second_y_status) {
  99.       if ((first_x_status == Positive_side) &&
  100.       (second_x_status == Positive_side))
  101.     num_of_crossings++;
  102.       else if ((first_x_status == Positive_side) || 
  103.            (second_x_status == Positive_side))
  104.     if ((u1 - v1*(u2-u1)/(v2-v1))>0)
  105.       num_of_crossings++;
  106.     }
  107.     first_x_status = second_x_status;    
  108.     first_y_status = second_y_status;    
  109.     u1 = u2; v1 = v2;
  110.   }
  111.   return(num_of_crossings&1); /* If odd -> is inside */
  112. }
  113.  
  114. /**********************************************************************/
  115. /* Form bounding box around a polygon                                 */
  116. /**********************************************************************/
  117. BoundingBoxType BoxPoly (poly)
  118.      Polygon *poly;
  119. {
  120.   BoundingBoxType box;
  121.   int i;
  122.  
  123.   box.min.x = box.max.x = poly->vert[0]->pos.x;
  124.   box.min.y = box.max.y = poly->vert[0]->pos.y;
  125.   box.min.z = box.max.z = poly->vert[0]->pos.z;
  126.   
  127.   for (i=1; i<poly->numVert; i++) {
  128.     if (poly->vert[i]->pos.x < box.min.x)
  129.       box.min.x = poly->vert[i]->pos.x;
  130.     else if (poly->vert[i]->pos.x > box.max.x)
  131.       box.max.x = poly->vert[i]->pos.x;
  132.     if (poly->vert[i]->pos.y < box.min.y)
  133.       box.min.y = poly->vert[i]->pos.y;
  134.     else if (poly->vert[i]->pos.y > box.max.y)
  135.       box.max.y = poly->vert[i]->pos.y;
  136.     if (poly->vert[i]->pos.z < box.min.z)
  137.       box.min.z = poly->vert[i]->pos.z;
  138.     else if (poly->vert[i]->pos.z > box.max.z)
  139.       box.max.z = poly->vert[i]->pos.z;
  140.   }
  141.   return (box);
  142. }
  143.  
  144. /**********************************************************************/
  145. /* Test if ray intersection just touches the poly edges               */
  146. /**********************************************************************/
  147. int RayPoly_Edge(r, pptr)
  148.      Vector *r;
  149.      Polygon *pptr;
  150. {
  151.   int i, ii;
  152.   int touch = FALSE;
  153.   Vector endpt[2];
  154.   Vector nrm;
  155.   static long PolyEdgeTests = 0;
  156.   static long PolyEdgeHits = 0;
  157.  
  158.   PolyEdgeTests++;
  159.   i= 0;
  160.   nrm = pptr->normal[0];
  161.   while((touch == FALSE) && (i < pptr->numVert)) {
  162.     ii = (i + 1) % pptr->numVert;
  163.     endpt[0] = pptr->vert[i]->pos;
  164.     endpt[1] = pptr->vert[ii]->pos;
  165.     if ((endpt[0].x <= r->x) && (r->x <= endpt[1].x))
  166.       if ((endpt[0].y <= r->y) && (r->y <= endpt[1].y))
  167.     if ((endpt[0].z <= r->z) && (r->z <= endpt[1].z))
  168.       if (dot(&nrm, vsub(r,&endpt[0])) == 0.0) {
  169.         touch = TRUE;
  170.         PolyEdgeHits++;
  171.       }
  172.     if ((endpt[0].x >= r->x) && (r->x >= endpt[1].x))
  173.       if ((endpt[0].y >= r->y) && (r->y >= endpt[1].y))
  174.     if ((endpt[0].z >= r->z) && (r->z >= endpt[1].z))
  175.       if (dot(&nrm, vsub(r,&endpt[0])) == 0.0) {
  176.         touch = TRUE;
  177.         PolyEdgeHits++;
  178.       }
  179.     i++;
  180.   }
  181.   return touch;
  182. }
  183.  
  184. /**********************************************************************/
  185. /* Ray plane intersect, return distance hit if new maximum */
  186. /**********************************************************************/
  187. int RayPlane(N,d, ray, hit)
  188.      Vector *N;
  189.      double d;
  190.      Ray *ray;
  191.      HitData *hit;
  192. {
  193.   double t;
  194.   static long PlaneTests = 0;
  195.   static long PlaneHits = 0;
  196.   Vector v;
  197.  
  198.   PlaneTests++;
  199.   
  200.   v = ray->dir;
  201.   norm(&v);
  202.   t = dot(N, &v);
  203.   if (fabs(t) < DAMN_SMALL)
  204.     return FALSE;
  205.   hit->distance = - (d + dot(N, &ray->origin)) / t;
  206.   hit->normal = *N;
  207.   PlaneHits++;
  208.   return TRUE;
  209. }
  210.  
  211. /**********************************************************************/
  212. /* Value of pt in plane equation */
  213. /**********************************************************************/
  214. #define Plane_Value(N,pt,d) (dot((N),(pt)) + (d))
  215.  
  216. /**********************************************************************/
  217. /* Test if point is within quadralateral or triangle */
  218. /**********************************************************************/
  219. int Inside_QuadTri(pptr, pt)
  220.      Polygon *pptr;
  221.      Vector *pt;
  222. {
  223.   int inside = FALSE;
  224.   polyUVW *uvw;
  225.   double pval;
  226.  
  227.   uvw = pptr->uvw;
  228.  
  229.   /* Test against U,V, W planes of traingle */
  230.   if (pptr->class == TRIANGLE) {
  231.     if (Plane_Value(&uvw->Nu, pt, uvw->du) >= 0.0)
  232.       if (Plane_Value(&uvw->Nv, pt, uvw->dv) >= 0.0)
  233.     if (Plane_Value(&uvw->Nw, pt, uvw->dw) >= 0)
  234.       inside = TRUE;
  235.   } 
  236.  
  237.   /* Test against U,V planes of quadralateral */
  238.   else {
  239.     pval = Plane_Value(&uvw->Nu, pt, uvw->du);
  240.     if (pval >= 0.0 && pval <= uvw->Lu) {
  241.       pval = Plane_Value(&uvw->Nv, pt, uvw->dv);
  242.       if (pval >= 0.0 && pval <= uvw->Lv) 
  243.     inside = TRUE;
  244.     }
  245.   }
  246.   return inside;
  247. }
  248.      
  249. /**********************************************************************/
  250. /* Ray polygon intersection. Return hit data, and if hit  */
  251. /**********************************************************************/
  252. int RayPolygon(pptr, ray, hit, max_distance, pbox, isect_quadtri)
  253.      Polygon *pptr;
  254.      Ray *ray;
  255.      HitData *hit;
  256.      double max_distance;
  257.      BoundingBoxType pbox;
  258.      int isect_quadtri;
  259. {
  260.   int hit_poly = FALSE;
  261.   double box_hit_dist;
  262.  
  263.   /* First check if box around polygon hits ray */
  264.   /* if (Bounds_Behind_Plane(&pbox, shootPatch->normal[0], shootPatch->d))
  265.      return FALSE;
  266.      */
  267.   if (!RayBoundingBox(pbox, *ray, &box_hit_dist))
  268.     return FALSE;
  269.   if (box_hit_dist >= max_distance)
  270.     return FALSE;
  271.  
  272.   /* Check plane to see if hits within a maximum distance.
  273.      If so, then check if inside polygon */
  274.   if (RayPlane(&pptr->normal[0], pptr->d, ray, hit)) {
  275.     if ((hit->distance > 0.0) && (hit->distance < max_distance)) {
  276.       Point_On_Ray(ray->origin,hit->distance,ray->dir,hit->intersect);
  277.       
  278.       if (isect_quadtri)          /* Intersect only quads and triangles */
  279.     hit_poly = Inside_QuadTri(pptr, &hit->intersect);
  280.       else {                      /* General poly intersect */
  281.     if ((hit_poly = RayPoly_Edge(&hit->intersect, pptr)) == FALSE)
  282.       hit_poly = Inside_Polygon(&hit->intersect, pptr);
  283.       }
  284.       /* Store hit info */
  285.       if (hit_poly) {
  286.     hit->obj = pptr->Mfather->Ofather;
  287.     hit->mesh = pptr->Mfather;
  288.     hit->poly = pptr;
  289.       }
  290.     }
  291.   }
  292.   return hit_poly;
  293. }
  294.  
  295. /**********************************************************************/
  296. /* Ray / polygonal mesh intersection. Return hit data, and if hit */
  297. /**********************************************************************/
  298. void RayMesh(ray, hit, optr)
  299.      register Ray *ray;
  300.      register HitData *hit;
  301.      register Objectt *optr;
  302. {
  303.   int i,j;
  304.   Mesh *mptr;
  305.   Polygon *pptr;
  306.   HitData polyhit;
  307.   BoundingBoxType pbox;
  308.  
  309.   /* Check all polygons of mesh with the ray */
  310.   mptr = optr->meshes;
  311.   i = 0;
  312.   while (i<optr->num_meshes && ray->visible > 0.0) {
  313.     RayStats.rayMesh++;
  314.     if (Option.debug)
  315.       printf("\t\tRay-Mesh: %s%d : visible = %g\n", 
  316.          mptr->name, mptr->id, ray->visible);
  317.     pptr = mptr->polys;
  318.     j = 0;
  319.     while (j<mptr->num_polys && ray->visible > 0.0) {
  320.       if (Option.debug)
  321.     printf("\t\tRay-Poly: %s%d : visible = %g\n", pptr->name, pptr->id,
  322.            ray->visible);
  323.  
  324.       /* Find if any ray / poly intersection */
  325.       polyhit = *hit;
  326.       if ((pptr != shootPatch) && (pptr != recPatch)) {
  327.     RayStats.rayPoly++;
  328.     pbox = BoxPoly(pptr);
  329.     if (RayPolygon(pptr, ray, &polyhit, hit->distance, pbox)) { 
  330.       if (polyhit.distance < hit->distance) { 
  331.         RayStats.intPoly++;
  332.         ray->visible = 0.0;
  333.         Store_HitInter(hit,polyhit.obj, polyhit.distance,
  334.                polyhit.intersect.x, polyhit.intersect.y,
  335.                polyhit.intersect.z);
  336.         if (!ray->shadow)
  337.           Store_HitNorm(hit, polyhit.normal.x, polyhit.normal.y,
  338.                 polyhit.normal.z);
  339.       }
  340.     }
  341.       }
  342.       pptr = pptr++; j++;
  343.     }
  344.     mptr = mptr++; i++;
  345.   }
  346. }
  347.